pacman::p_load(sf, spdep, sfdep, tmap, tidyverse, plotly, Kendall)Take-home Exercise 3: Quake Quest - EDA + Spatial Autocorrelation of
Project Motivation
Given Indonesia’s geographical location along the Pacific Ring of Fire, the country experiences frequent and intense earthquakes. It lies in the intersections with three major tectonic plates: the Indo-Australian Plate, the Eurasian Plate, and the Pacific Plate. The constant interaction of these plates, along with the interaction with several other minor plates such as the Philippine Sea plate and the Caroline plate, amplifies the geological complexity of the region. This makes it crucial to examine the geospatial data related to seismic activity. We aim to also use this data to conduct further research to understand how geological features might relate to these events.
Earthquakes in Indonesia often strike in densely populated areas, causing many people to lose their homes, livelihoods and property. By leveraging spatial analysis techniques, we can provide insights that can inform proactive disaster management strategies, urban planning decisions and infrastructure resilience in Indonesia. Keeping well-informed is a good way to minimise negative impacts in the event of a seismic activity.
Furthermore, earthquakes are often met with foreshocks and aftershocks of smaller magnitudes, which can be analysed to understand the spatial patterns and predict coming shocks. This includes identifying areas prone to high levels of seismic activity and implementing building codes and infrastructure improvements to reduce the impact of earthquakes on communities.
Performing such spatial analysis necessitates expertise in execution, including proficiency in data wrangling and cleaning to ensure accuracy. However, many users lack coding skills or familiarity with R packages, posing a challenge to understanding and conducting these analyses effectively. Consequently, developing an interactive application offers an accessible solution. Users can comprehend analysis findings and adjust parameters to suit their objectives without the lengthy execution processes.
Project Objectives
Conducting Exploratory Data Analysis (EDA) provides users with contextual insights into the regions where frequent seismic activities occur.
Conduct spatial autocorrelation analysis
Global and Local Moran’s I statistical test to examine the degree of similarity between observations in a dataset based on their spatial proximity
Identify specific areas of significant spatial clustering (hotspots and cold spots)
- Spatio-temporal analysis to identify emerging hotspots
Data Source
This Earthquake dataset is taken from the Earthquake Repository managed by BMKG (an Indonesian non-departmental government agency). It contains earthquake event data from 1 Nov 2008 to 15 Dec 2023, but may not be accurate for some of the last earthquake events recorded.
Indonesia’s Subnational Administrative Boundaries extracted from HDX.
Checking data
Importing Geospatial Data
indonesia <- st_read(dsn = "data/geospatial/idn_adm_bps_20200401_shp",
layer = "idn_admbnda_adm3_bps_20200401")Reading layer `idn_admbnda_adm3_bps_20200401' from data source
`/Users/shanellelam/Library/CloudStorage/OneDrive-SingaporeManagementUniversity/School/Year 3/chanelelele/IS415-GAA/Take-home_Ex/Take-home_Ex03/data/geospatial/idn_adm_bps_20200401_shp'
using driver `ESRI Shapefile'
Simple feature collection with 7069 features and 16 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 95.01079 ymin: -11.00762 xmax: 141.0194 ymax: 6.07693
Geodetic CRS: WGS 84Extracting only West Java regions
west_java <- indonesia[indonesia$ADM1_EN == "Jawa Barat", ]west_java <- st_transform(west_java, crs = 23830)plot(st_geometry(west_java))
indo_earthquake <- read_csv("data/aspatial/katalog_gempa.csv")head(indo_earthquake)# A tibble: 6 × 13
tgl ot lat lon depth mag remark strike1 dip1 rake1 strike2
<date> <time> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <dbl>
1 2008-11-01 21:02:43 -9.18 119. 10 4.9 Sumba… NA NA NA NA
2 2008-11-01 20:58:50 -6.55 130. 10 4.6 Banda… NA NA NA NA
3 2008-11-01 17:43:12 -7.01 107. 121 3.7 Java … NA NA NA NA
4 2008-11-01 16:24:14 -3.3 128. 10 3.2 Seram… NA NA NA NA
5 2008-11-01 16:20:37 -6.41 130. 70 4.3 Banda… NA NA NA NA
6 2008-11-01 14:47:00 -7.37 105. 18 3.3 Java … NA NA NA NA
# ℹ 2 more variables: dip2 <dbl>, rake2 <dbl># Convert to sf object
indoEarthq_sf <- st_as_sf(indo_earthquake, coords = c("lon", "lat"), crs = "+proj=longlat +datum=WGS84") # Transform the geometry to EPSG:23830
indoEarthq_sf <- st_transform(indoEarthq_sf, crs = "+init=EPSG:23830")indoEarthq_filter <- indoEarthq_sf %>%
filter(tgl >= as.Date("2009/01/01") & tgl <= as.Date("2023/12/31"))# Group by the 'remark' column and calculate the count of each type
remark_counts <- indoEarthq_filter %>% group_by(remark) %>%
summarise(count = n())
# View the resulting count of each type in the 'remark' column # Sort the data frame by count in descending order and select the top 10 rows
top_10_remark <- remark_counts %>%
arrange(desc(count)) %>%
head(10)
# Plot a bar graph
ggplot(top_10_remark,
aes(x = remark, y = count)) +
geom_bar(stat = "identity", fill = "skyblue") +
labs(title = "Top 10 Region Eathquake Count",
x = "Region",
y = "Count") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot(st_geometry(indoEarthq_sf))
java_earthq <- st_intersection(indoEarthq_filter, west_java)
java_earthqSimple feature collection with 1382 features and 27 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 1522045 ymin: 614269.8 xmax: 1792451 ymax: 816381.6
Projected CRS: DGN95 / Indonesia TM-3 zone 46.2
# A tibble: 1,382 × 28
tgl ot depth mag remark strike1 dip1 rake1 strike2 dip2
* <date> <time> <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 2015-03-13 17:23:27 15 2.9 Java - Ind… NA NA NA NA NA
2 2016-01-11 16:58:36 18 2.9 Java - Ind… NA NA NA NA NA
3 2016-11-30 08:00:24 60 3.2 Java - Ind… NA NA NA NA NA
4 2017-08-09 04:57:31 111 3.6 Java - Ind… NA NA NA NA NA
5 2018-01-03 03:49:07 110 3.3 Java - Ind… NA NA NA NA NA
6 2018-01-26 12:40:02 96 3 Java - Ind… NA NA NA NA NA
7 2020-03-22 20:27:05 91 3.5 Java - Ind… NA NA NA NA NA
8 2022-09-03 02:29:58 109 2.5 Java - Ind… NA NA NA NA NA
9 2023-01-11 05:20:19 16 2.5 Java - Ind… NA NA NA NA NA
10 2015-03-28 10:25:24 10 3.3 Java - Ind… NA NA NA NA NA
# ℹ 1,372 more rows
# ℹ 18 more variables: rake2 <dbl>, Shape_Leng <dbl>, Shape_Area <dbl>,
# ADM3_EN <chr>, ADM3_PCODE <chr>, ADM3_REF <chr>, ADM3ALT1EN <chr>,
# ADM3ALT2EN <chr>, ADM2_EN <chr>, ADM2_PCODE <chr>, ADM1_EN <chr>,
# ADM1_PCODE <chr>, ADM0_EN <chr>, ADM0_PCODE <chr>, date <date>,
# validOn <date>, validTo <date>, geometry <POINT [m]>Point with Polygon Join: If west_java contains polygon geometries and java_earthq contains point geometries, and you perform a spatial join to find points within polygons, the resulting dataset will likely have polygon geometries. This is because the join operation will match points to polygons, resulting in the polygons being retained.
java_eq_left <- st_join(west_java, java_earthq) %>%
rename(ADM3_EN = ADM3_EN.x)unique_values <- unique(java_eq_left$ADM2_EN.x)
print(unique_values) [1] "Cianjur" "Kota Bandung" "Indramayu" "Majalengka"
[5] "Bandung" "Cirebon" "Bogor" "Purwakarta"
[9] "Bekasi" "Kota Banjar" "Ciamis" "Garut"
[13] "Sukabumi" "Kota Bekasi" "Tasikmalaya" "Karawang"
[17] "Kota Sukabumi" "Bandung Barat" "Kota Depok" "Subang"
[21] "Kota Bogor" "Sumedang" "Kota Tasikmalaya" "Kuningan"
[25] "Pangandaran" "Kota Cimahi" "Kota Cirebon" "Waduk Cirata" tmap_mode('view')
# Define basemaps
osm_basemap <- tm_basemap(server = "OpenStreetMap.HOT")
imagery_basemap <- tm_basemap(server = "Esri.WorldImagery")
# Plot dl_class over OpenStreetMap basemap
eo_map <- osm_basemap +
tm_shape(java_eq_left) +
tm_markers(col = "red", size = 1) + # Adjust marker color and size as needed
tmap_options(check.and.fix = TRUE)
# View the map
eo_mapFor the following EDA map, the user will be able to tailor their exploration by filtering seismic events based on specific geographical regions – provinces and city regions.
Paramters to utilize the EDA map:
- Province Selection:
- Users can input the name of the province they wish to investigate into the designated text input field.
- This province selection aligns with the geographical delineation provided in the “ADM2_EN” column of our dataset.
- City Region Specification:
- For a more granular analysis, users can further refine their search by inputting the name of the city region they are interested in.
- Date Range:
- In addition, users have the flexibility to specify the date range of their interest, allowing them to narrow down their analysis to seismic events occurring within a particular timeframe
- The city region input corresponds to the localized districts recorded in the “ADM3_EN” column of our dataset.
Navigating Filtered Earthquake Data:
Upon entering the desired province and city region, the map dynamically updates to exclusively display seismic events recorded within the defined area and timeframe. Each marker on the map signifies the precise location of an earthquake event, providing users with a visual representation of seismic activity within their specified regions of interest.
tmap_mode('plot')
tm_shape(west_java)+
tm_polygons()+
tm_shape(java_earthq)+
tm_dots()
java_earthq_count <- java_eq_left %>%
group_by(ADM3_EN) %>%
summarize(num_rows_in_group = n())# Plot choropleth map
tm_shape(java_earthq_count) +
tm_polygons("num_rows_in_group", palette = "Blues", title = "Earthquake Count") +
tm_layout(title = "Choropleth Map of Earthquake Count in West Java")
Global Spatial Autocorrelation
Spatial autocorrelation is the term used to describe the presence of systematic spatial variation in a variable.
Hypotheses:
Null Hypothesis: The map suggests that the seismic activity in West Java occurred randomly over space. (complete spatial randomness)
Alternative Hypothesis: The map suggests that the seismic activity in West Java occurred in a non-random pattern over space (spatially clustered)
wm_q <- java_earthq_count %>%
mutate(nb = st_contiguity(geometry),
wt = st_weights(nb,
style = "W"),
.before = 1) The st_weights() function in the sf package is used to compute spatial weights for spatial analysis. It generates weights for spatial entities based on their relationships with neighboring entities.
Paramters to generate the LISA map: Here’s an overview of the options available for the style parameter in st_weights():
W (Row Standardized Weights): This is the default option. It computes weights that are row standardized, meaning the sum of weights over all links connected to a particular spatial entity sums to the entity’s row number (n).
B (Basic Binary Coding): This option generates basic binary weights where neighbors are assigned a weight of 1 if they share a boundary with the entity, and 0 otherwise.
C (Globally Standardized Weights): It computes weights that are globally standardized, meaning the sum of weights over all links sums to the total number of neighbors (n).
U (Equal to C Divided by the Number of Neighbors): This option is similar to C, but the weights are divided by the number of neighbors, resulting in a sum of weights over all links equal to unity. -
minmax: This option computes weights such that the maximum weight of any neighbor is set to 1, and all other weights are scaled accordingly.
S (Variance-Stabilizing Coding Scheme): This option applies a variance-stabilizing transformation to the weights, which can help stabilize variance and mitigate the impact of outliers in spatial analysis.
wm_qSimple feature collection with 584 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 1519914 ymin: 609162.8 xmax: 1802448 ymax: 829713.7
Projected CRS: DGN95 / Indonesia TM-3 zone 46.2
# A tibble: 584 × 5
nb wt ADM3_EN num_rows_in_group geometry
* <nb> <list> <chr> <int> <MULTIPOLYGON [m]>
1 <int [6]> <dbl [6]> Agrabinta 9 (((1585322 671054.1, 1585…
2 <int [9]> <dbl [9]> Andir 1 (((1656295 717797.2, 1656…
3 <int [5]> <dbl [5]> Anjatan 1 (((1702273 781735.2, 1702…
4 <int [4]> <dbl [4]> Antapani 1 (((1667227 716133.9, 1667…
5 <int [4]> <dbl [4]> Arahan 1 (((1732515 778991.6, 1732…
6 <int [7]> <dbl [7]> Arcamanik 1 (((1667805 716062.5, 1667…
7 <int [9]> <dbl [9]> Argapura 1 (((1742445 718654.8, 1742…
8 <int [4]> <dbl [4]> Arjasari 2 (((1662207 702189.3, 1662…
9 <int [8]> <dbl [8]> Arjawinangun 1 (((1750909 748231.8, 1750…
10 <int [6]> <dbl [6]> Astanaanyar 1 (((1660461 714396.9, 1660…
# ℹ 574 more rowsComputing Global Moran I
moranI <- global_moran(wm_q$num_rows_in_group,
wm_q$nb,
wm_q$wt)
glimpse(moranI)List of 2
$ I: num 0.301
$ K: num 244Performing Global Moran’sI test
global_moran_test(wm_q$num_rows_in_group,
wm_q$nb,
wm_q$wt)
Moran I test under randomisation
data: x
weights: listw
Moran I statistic standard deviate = 16.297, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.3006401204 -0.0017152659 0.0003442278 Performing Global Moran’I permutation test
In practice, monte carlo simulation should be used to perform the statistical test. For sfdep, it is supported by globel_moran_perm() It is alway a good practice to use set.seed() before performing simulation. This is to ensure that the computation is reproducible.
set.seed(1234)global_moran_perm(wm_q$num_rows_in_group,
wm_q$nb,
wm_q$wt,
nsim = 99)
Monte-Carlo simulation of Moran I
data: x
weights: listw
number of simulations + 1: 100
statistic = 0.30064, observed rank = 100, p-value < 2.2e-16
alternative hypothesis: two.sidedComputing local Moran’s I
lisa <- wm_q %>%
mutate(local_moran = local_moran(
num_rows_in_group, nb, wt, nsim = 99),
.before = 1) %>%
unnest(local_moran)Visualising local Moran’s I
tmap_mode("plot")
tm_shape(lisa) +
tm_fill("ii") +
tm_borders(alpha = 0.5) +
tm_view(set.zoom.limits = c(6,8)) +
tm_layout(main.title = "local Moran's I of sesismic activity",
main.title.size = 0.8)
Visualising p-value of local Moran’s I
tmap_mode("plot")
tm_shape(lisa) +
tm_fill("p_ii_sim") +
tm_borders(alpha = 0.5) +
tm_layout(main.title = "p-value of local Moran's I",
main.title.size = 0.8)
Visuaising local Moran’s I and p-value
Added the variable midpoint=NA to map1 as the variable(s) “ii” contains positive and negative values, which causes midpoint to be set to 0. Therefore, I set midpoint to NA to show the full spectrum of the color palette.
tmap_mode("plot")
map1 <- tm_shape(lisa) +
tm_fill("ii", midpoint = NA) +
tm_borders(alpha = 0.5) +
tm_view(set.zoom.limits = c(6,8)) +
tm_layout(main.title = "local Moran's I of sesismic activity",
main.title.size = 0.8)
map2 <- tm_shape(lisa) +
tm_fill("p_ii",
breaks = c(0, 0.001, 0.01, 0.05, 1),
labels = c("0.001", "0.01", "0.05", "Not sig")) +
tm_borders(alpha = 0.5) +
tm_layout(main.title = "p-value of local Moran's I",
main.title.size = 0.8)
tmap_arrange(map1, map2, ncol = 2)
Renaming the column from ADM3_EN to City Region:
lisa <- lisa %>%
rename(`City Region` = ADM3_EN)## added the tooltip
tmap_mode("plot")
map1 <- tm_shape(lisa) +
tm_fill("ii", popup.vars = c("City Region"), midpoint = NA) +
tm_borders(alpha = 0.5) +
tm_view(set.zoom.limits = c(6,8)) +
tm_layout(main.title = "Local Moran's I of seismic activity",
main.title.size = 0.8)
# Map 2: p-value of local Moran's I
map2 <- tm_shape(lisa) +
tm_fill("p_ii",
popup.vars = c("City Region"),
breaks = c(0, 0.001, 0.01, 0.05, 1),
labels = c("0.001", "0.01", "0.05", "Not sig")) +
tm_borders(alpha = 0.5) +
tm_layout(main.title = "p-value of Local Moran's I",
main.title.size = 0.8)
tmap_arrange(map1, map2, ncol = 2)
The added tooltip allows the user to hover over the sector and check the name of the region.
Visualising LISA map
# attempt 1
lisa_sig <- lisa %>%
filter(p_ii < 0.05)
tmap_mode("plot")
tm_shape(lisa) +
tm_polygons() +
tm_borders(alpha = 0.5) +
tm_shape(lisa_sig) +
tm_fill("mean") +
tm_borders(alpha = 0.4)
tmap_mode("plot")
map <- tm_shape(lisa) +
tm_polygons() +
tm_borders(alpha = 0.5) +
tm_shape(lisa_sig) +
tm_fill("mean", popup.vars = "City Region") +
tm_borders(alpha = 0.4)
map
The added tooltip allows the user to hover over the sector and check the name of the region.
The LISA map is an interpreted map by combining local Moran’s I of geographical areas and their respective p-values.
Paramters to generate the LISA map:
The fields contain the LISA categories:
mean: Represents the average value of the variable of interest (e.g., seismic activity) within the spatial neighborhood of each entity. It is computed as the mean value of the variable across neighboring entities.
median: Represents the median value of the variable within the spatial neighborhood of each entity. Unlike the mean, which is influenced by extreme values or outliers, the median provides a measure of central tendency that is less sensitive to outliers.
PySAL: Represents the local spatial autocorrelation statistic calculated using the PySAL library, which is a Python library for spatial analysis. This statistic measures the degree of spatial clustering or dispersion of the variable of interest within the spatial neighborhood of each entity. It indicates whether the local pattern is similar to the global pattern (positive spatial autocorrelation), dissimilar (negative spatial autocorrelation), or random (no spatial autocorrelation).
Hot Spot and Cold Spot Area Analysis (HCSA)
HCSA uses spatial weights to identify locations of statistically significant hot spots and cold spots in an spatially weighted attribute that are in proximity to one another based on a calculated distance. The analysis groups features when similar high (hot) or low (cold) values are found in a cluster. The polygon features usually represent administration boundaries or a custom grid structure.
Computing local Gi* statistics
# Assuming java_eq_left is an sf object with geometry column named "geometry"
java_earthq_count <- java_earthq_count %>%
mutate(geometry_point = st_centroid(geometry)) # Convert polygon geometries to point geometries (centroids)wm_idw <- java_earthq_count %>%
mutate(nb = st_contiguity(geometry),
wts = st_inverse_distance(nb, geometry_point,
scale = 1,
alpha = 1),
.before = 1)Parameters for Hot Spot and Cold Spot Area Analysis (HCSA):
- Scale:
- Scale parameter in the st_inverse_distance function that controls the scaling factor for the inverse distance weights.
- Adjusting this parameter will impact the influence of distant neighbors on the computed weights.
- Higher values of scale result in a more gradual decrease in weight with distance, while lower values give more weight to closer neighbors.
- Alpha:
- The alpha parameter in the st_inverse_distance function controls how fast the weight decreases as distance increases.
- A higher alpha value makes the weight decrease faster, while lower values make it decrease slower.
- Number of Simulations (nsim):
- In the local_gstar_perm function, the nsim parameter specifies the number of permutations to use for the Monte Carlo simulation.
- Increasing nsim can lead to more accurate p-values. However, it also increases computation time.
By adjusting these parameters, users can tailor their exploration of hot spot and cold spot areas in the seismic data, allowing for a more nuanced analysis of spatial patterns and relationships.
HCSA <- wm_idw %>%
mutate(local_Gi = local_gstar_perm(
num_rows_in_group, nb, wt, nsim = 99),
.before = 1) %>%
unnest(local_Gi)
HCSASimple feature collection with 584 features and 12 fields
Active geometry column: geometry
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 1519914 ymin: 609162.8 xmax: 1802448 ymax: 829713.7
Projected CRS: DGN95 / Indonesia TM-3 zone 46.2
# A tibble: 584 × 14
gi_star e_gi var_gi p_value p_sim p_folded_sim skewness kurtosis nb
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <nb>
1 1.14 0.00246 0.00000639 0.628 0.530 0.24 0.12 3.86 <int>
2 -0.671 0.00170 0.00000280 -0.678 0.498 0.08 0.01 4.13 <int>
3 -0.518 0.00174 0.00000471 -0.539 0.590 0.16 0.01 5.09 <int>
4 -0.472 0.00182 0.00000718 -0.468 0.640 0.3 0.01 5.01 <int>
5 -0.426 0.00131 0.00000142 -0.533 0.594 0.78 0.39 3.12 <int>
6 -0.599 0.00141 0.00000157 -0.671 0.502 0.16 0.01 3.00 <int>
7 -0.306 0.00157 0.00000215 -0.260 0.795 0.92 0.47 3.87 <int>
8 -0.192 0.00160 0.00000543 -0.152 0.879 0.74 0.41 7.03 <int>
9 -0.636 0.00148 0.00000150 -0.744 0.457 0.08 0.01 2.61 <int>
10 -0.560 0.00189 0.00000777 -0.475 0.635 0.12 0.01 4.43 <int>
# ℹ 574 more rows
# ℹ 5 more variables: wts <list>, ADM3_EN <chr>, num_rows_in_group <int>,
# geometry <MULTIPOLYGON [m]>, geometry_point <POINT [m]>tmap_mode("plot")
tm_shape(HCSA) +
tm_fill("gi_star") +
tm_borders(alpha = 0.5) +
tm_view(set.zoom.limits = c(6,8))
tmap_mode("plot")
tm_shape(HCSA) +
tm_fill("p_sim") +
tm_borders(alpha = 0.5)
tmap_mode("plot")
map1 <- tm_shape(HCSA) +
tm_fill("gi_star") +
tm_borders(alpha = 0.5) +
tm_view(set.zoom.limits = c(6,8)) +
tm_layout(main.title = "Gi* of GDPPC",
main.title.size = 0.8)
map2 <- tm_shape(HCSA) +
tm_fill("p_value",
breaks = c(0, 0.001, 0.01, 0.05, 1),
labels = c("0.001", "0.01", "0.05", "Not sig")) +
tm_borders(alpha = 0.5) +
tm_layout(main.title = "p-value of Gi*",
main.title.size = 0.8)
tmap_arrange(map1, map2, ncol = 2)
Visualising hot spot and cold spot areas
HCSA_sig <- HCSA %>%
filter(p_sim < 0.05)
tmap_mode("plot")
map <- tm_shape(HCSA) +
tm_polygons() +
tm_borders(alpha = 0.5) +
tm_shape(HCSA_sig) +
tm_fill("gi_star", popup.vars = "ADM3_EN") +
tm_borders(alpha = 0.4)
map
The added tooltip allows the user to hover over the sector and check the name of the region.
class(HCSA)[1] "sf" "tbl_df" "tbl" "data.frame"Storyboard
